{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qmt.physics_constants import parse_unit,to_float\n",
    "from qmt.geometry import part_3d, build_3d_geometry\n",
    "from qmt.materials import Materials, build_materials, make_materials_library\n",
    "import pickle, logging\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.interpolate import griddata\n",
    "import deepdish\n",
    "from qms.physics import PoissonFem\n",
    "import fenics as fn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load some 2D shapes from simple_wire.FCStd, and then build some a 3D structure out of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "back_gate = part_3d.ExtrudePart(\"back_gate\", \"Sketch\", z0=-0.2, thickness=0.1)\n",
    "vacuum = part_3d.ExtrudePart(\"vacuum\", \"Sketch003\", z0=-0.5, thickness=1.0)\n",
    "wire = part_3d.ExtrudePart(\"wire\", \"Sketch002\", z0=0.0, thickness=0.1)\n",
    "shell = part_3d.ExtrudePart(\"shell\", \"Sketch002\", z0=0.1, thickness=0.05)\n",
    "build_order = [wire, shell, back_gate, vacuum]\n",
    "file_path = './simple_wire.FCStd'\n",
    "geo_data = build_3d_geometry(input_parts=build_order, input_file=file_path,\n",
    "                             xsec_dict={'central':{'axis':(1.,0.,0.),'distance':0.}})\n",
    "geo_data.write_fcstd('built_geo.fcstd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point you can try opening built_geo.fcstd with FreeCAD and taking a look at the built shape. Feel free to skip this step if you're unfamiliar with the FreeCAD GUI.\n",
    "We can check that our cross sections and parts are as expected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_data.xsecs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geo_data.parts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can do meshing. A mesh has already been pre built so you can just load it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qms.meshing import MeshPart, MeshData\n",
    "# from qms.tasks import build_3d_mesh\n",
    "# mesh_data = build_3d_mesh(\n",
    "#     geo_data,\n",
    "#     {\n",
    "#         \"back_gate\": MeshPart(mesh_max_size=0.1),\n",
    "#         \"vacuum\": MeshPart(mesh_max_size=0.05),\n",
    "#         \"wire\": MeshPart(mesh_max_size=0.01),\n",
    "#         \"shell\": MeshPart(mesh_max_size=0.01),\n",
    "#     },\n",
    "#     \"comsol\"\n",
    "# )\n",
    "\n",
    "# mesh_data.save(\"mesh_data.h5\")\n",
    "\n",
    "mesh_data = MeshData.load(\"mesh_data.h5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the reference level and Al work function to get a 0.1 meV band offset between InSb and Al:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_lib = Materials()\n",
    "Al_WF = mat_lib['Al']['workFunction']\n",
    "InSb_EA = mat_lib['InSb']['electronAffinity']\n",
    "InSb_BG = mat_lib['InSb']['directBandGap']\n",
    "InSb_VBO = mat_lib['InSb']['valenceBandOffset']\n",
    "Al_WF_level = 0.0-(Al_WF)\n",
    "InSb_CB_level = 0.0-InSb_EA+InSb_VBO\n",
    "WF_shift = 200.*parse_unit('meV')-(Al_WF_level-InSb_CB_level)\n",
    "new_Al_WF = (Al_WF-WF_shift)\n",
    "ref_level = -new_Al_WF\n",
    "\n",
    "mat_lib = make_materials_library({\"Al\":{\"workFunction\": new_Al_WF}})\n",
    "mat_data = build_materials(geo_data, \n",
    "                           {\"back_gate\": \"Al\", \"vacuum\": \"air\", \"wire\": \"InSb\", \"shell\": \"Al\"},\n",
    "                           mat_lib)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the electrostatic simulations. Again you can just load the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from qms.tasks.thomas_fermi import ThomasFermiPart\n",
    "# from qms.tasks import run_3d_thomas_fermi\n",
    "# import sympy.physics.units as spu\n",
    "# logging.basicConfig(level=logging.INFO)\n",
    "# tf_data = run_3d_thomas_fermi(geo_data,\n",
    "#                               mesh_data,\n",
    "#                               mat_data,\n",
    "#                               {\n",
    "#                                   \"back_gate\": ThomasFermiPart(\"metal_gate\", boundary_condition={\"voltage\": 0.0 * spu.V}),\n",
    "#                                   \"vacuum\": ThomasFermiPart(\"dielectric\"),\n",
    "#                                   \"wire\": ThomasFermiPart(\"semiconductor\"),\n",
    "#                                   \"shell\": ThomasFermiPart(\"metal_gate\", boundary_condition={\"voltage\": 0.0 * spu.V}),\n",
    "#                               },\n",
    "#                               reference_level=ref_level,\n",
    "#                               order=1,\n",
    "#                              )\n",
    "# tf_data.save(\"tf_data.h5\")\n",
    "\n",
    "from qms.physics import PoissonFem\n",
    "tf_data = PoissonFem.load(\"tf_data.h5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Let's take a look at the potential profile to make sure it looks reasonable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "femcoords = tf_data.coordinates\n",
    "femdata = tf_data.potential\n",
    "xgrid, ygrid, zgrid = np.mgrid[0:1, -0.2:0.2:0.005, -0.5:0.5:0.0125]\n",
    "plot_potential = griddata(femcoords, femdata, (xgrid, ygrid, zgrid), method='linear')\n",
    "plt.pcolor(ygrid[0],zgrid[0],plot_potential[0])\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks fine. Let's now look at a line cut:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgrid, ygrid, zgrid = np.mgrid[0:1, 0:1, -0.2:0.2:0.002]\n",
    "potential_cut = griddata(femcoords, femdata, (xgrid, ygrid, zgrid), method='linear')\n",
    "plt.plot(zgrid[0,0],potential_cut[0,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This all looks fine. In the region between 0.0 and 0.1, we have accumulation. Let's make sure this holds up when taking into account the conduction band offset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgrid, ygrid, zgrid = np.mgrid[0:1, 0:1, 0:0.1:0.0005]\n",
    "potential_cut = griddata(femcoords, femdata, (xgrid, ygrid, zgrid), method='linear')\n",
    "plt.plot(zgrid[0,0],potential_cut[0,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting in the proper band offsets, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zvec = zgrid[0,0]\n",
    "potential_cut = potential_cut[0,0]\n",
    "offset_CB = to_float((InSb_CB_level-ref_level) /parse_unit('meV'))/1e3\n",
    "offset_VB = offset_CB-InSb_BG/parse_unit('meV')/1e3\n",
    "plt.plot(zvec,offset_CB-potential_cut)\n",
    "plt.plot(zvec,offset_VB-potential_cut)\n",
    "plt.plot(zvec,np.zeros(zvec.shape))\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
